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projector method to determine the topological susceptibility in the pure Yang-Mills theory. 
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mass fermions including for the hrst time also the flavour singlet decay constant. The 
Witten-Veneziano formula represents a leading order expression in the framework of chiral 
perturbation theory and we also employ leading order chiral perturbation theory to relate 
the flavor singlet decay constant to the relevant decay constant parameters in the quark 
flavor basis and flavor non-singlet decay constants. After taking the continuum and the 
SU(2) chiral limits we compare both sides and hnd good agreement within uncertainties. 
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1 Introduction 


The Witten-Veneziano formula [1, 2] was first derived by Witten for massless quarks 

V," = ■ (I'l) 

h 

It aims at providing an explanation for the unexpectedly large mass of the t]' meson, by 
relating its mass in the chiral limit M^i to non-trivial topological fluctuations of the gauge 
fields, which are encoded in the topological susceptibility computed in pure Yang-Mills 
(YM) theory Xoo- Furthermore, the formula involves the singlet decay constant /o and the 
number of quark flavors Nf. 

In order to obtain Eq. (1.1) one considers the limit of a large number of colors Nc, for 
which simplifications of the theory occur that allow one to address a variety of problems 
which are otherwise impossible to investigate. In particular, in the’t Hooft limit (Nc —)• oo, 
while g'^Nc and Nf are kept fixed, with g denoting the gauge coupling) the anomalously 
broken axial 17(1) symmetry is restored and the rj' becomes a Goldstone boson, i.e. Mj^i —)• 0 
[1, 3]. The formula itself is valid up to corrections of 0{1/Nc) while and l//o vanish 
as 0{1/Nc). However, the topological susceptibility which appears on the right-hand side 
of Eqs. (1.1) remains finite in the large W limit, i.e. Xoo = 0(1)) such that in the W oo 
limit both sides of Eq. (1.1) vanish. The question remains whether Nc = 3 in QCD is large 
enough in practice to sufficiently suppress corrections to this formula to higher order in 
1/Nc. 

An alternative way to derive the Witten-Veneziano formula is to expand the anomalous 
flavor-singlet Ward-Takahashi identities of the theory order by order in n = Nf /Nc around 
n = 0 [2, 4]. Again, the Witten-Veneziano formula corresponds to the lowest order relation 
in this expansion. Besides, we remark that for lattice QCD it is also possible to obtain an 
unambiguous, theoretical sound implementation of the Witten-Veneziano formula through 
the study of anomalous flavor-singlet Ward-Takahashi identities in the limit rt —>■ 0 [5, 6]. 

However, the very first attempts to compute that topological susceptibility in YM 
theory date back to the early 80s [7, 8] and a first reasonable number was found in [9]. 

Moving away from the chiral limit leads to corrections to the formula which are linear 
in the quadratic meson masses 

^(M2 + m2 - 2Ml) = Xoo . (1.2) 

This result was first derived in [2] from the aforementioned expansion in u. In the above 
expression terms have been reshuffled compared to Eq. (1.1) to isolate the topological 
susceptibility in the YM theory on the r.h.s., whereas the l.h.s. of the formula contains 
only quantities that need to be computed in full QCD. 

From a modern point of view, the Witten-Veneziano formula can also be derived from 
effective field theory employing a combined power counting scheme in quark masses niq, 
momenta p and l/A^o given by niq = 0{5), p^O{5) and l/W = 0(6), where <5 is a small 
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expansion parameter [10-13]. In principle, this approach allows one to systematically cal¬ 
culate higher order corrections to the standard form of Eq. (1.2) 

9 9 9 AN f 

Mrj -I- M^j — 2 M^ = —^Xoo , (1-3) 

JO 

which represents the leading order expression with respect to the expansion in 6 used in 
chiral perturbation theory (xPT). Note that in general /o 7 ^ f-rr holds, as the singlet decay 
constant /o in Eqs. (1.1), (1-2) becomes equal to the octet decay constant /g only when 
simultaneously taking the chiral limit and dropping all corrections in 1/Nc [11]. Since on 
the lattice we do not work in the octet-singlet basis, we cannot compute /o directly but need 
to apply leading order chiral perturbation theory in order to relate fo to the corresponding 
decay constant parameter in the quark flavor basis and the relevant flavor non-singlet decay 
constants /jr and fjc- 

In this paper, we compute both the topological susceptibility in the pure YM theory (or 
quenched QCD) and the meson masses and the singlet decay constant fo in full QCD. After 
discussing the lattice actions, we will first discuss the determination of Xoo using the so- 
called spectral projector method [14] based on dedicated quenched simulations^. Thereafter, 
we discuss the determination of M^, M^/, and in particular fo using Nf = 2 -|- 1 -|- 1 
lattice QCD. Finally, we compare both results for Xoo, finding good agreement between 
quenched and dynamic computation. 


^For early attempts to compute the topological susceptibility in YM theory we refer to [7, 8], while a 
first reasonable number was found in [9]. 
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2 Lattice Actions 


2.1 Nf = 2 + 1 + 1 lattice QCD 

The calculation of the masses presented in this work was performed using configurations 
with Nf = 2 + 1 + 1 dynamical Wilson twisted mass fermions at maximal twist generated 
by the European Twisted Mass Collaboration (ETMC) [15-17]. 

The lattice Wilson twisted mass fermion action^ for the light sector [19-21], i.e. u and 
d quarks, is given, in the twisted mass basis, by: 

= a^'^Xi{x){Dw + mo + ii2a5T3)xi{x ), ( 2 . 1 ) 

X 

where mo and m are the bare untwisted and twisted quark masses, respectively. Xl = 
iXu Xd)'^ is a flavor doublet and ra acts in flavor space. The massless Wilson Dirac operator 
Dw is defined as 

^w^ = ^(7M(V^ + V*)-aV*V^), (2.2) 

where and V* are the forward and backward covariant derivatives. 

The lattice Wilson twisted mass action for the heavy doublet Xh = (XcXs)^ [21, 22], 
i.e. strange and charm quarks, is given by: 


Sh[il^,ip,U] = a^'^Xh{x){Dw + mo + + ^i5n)Xh{x) ■ (2.3) 

X 


The bare twisted mass parameters and /ig- are related to the bare strange and charm 
quark masses through the following relation: 


+c,s — 


d-a i +5 


(2.4) 


where Zp jZs dehnes the ratio of pseudoscalar and scalar flavor non-singlet renormalization 
factors. Both doublets, xi and Xh, are related to their counterparts in the physical basis 
via chiral rotations. 

For the gauge sector the Iwasaki [23, 24] gauge action was used which is defined as: 


Sg[U] = ^0 E I^^Tr (1 - + 6i 


/ 


V 




4 

E 


ReTV (1 - P^: 


lx2\ 


(2.5) 


where bi = —0.331 and 6o = 1 ~ 86i. The dynamical ensembles used in this work are 
compiled in Tab. 1 with the labeling adopted from Ref. [15]. For each /3-value, there 
are several quark masses available, which allows one to study the chiral extrapolation. 
The input parameters of our simulations are supplemented by Tab. 2, which contains the 
chirally extrapolated values of the Sommer parameter ro/a for each value of /3, which we 
use to set the scale in our dynamical simulations. In addition, we quote in this table the 
chirally extrapolated values of the renormalization constant ratio Z = ZpjZs from the two 
dehnitions discussed in Ref. [25]. These two definitions differ by lattice artefacts and are, 
therefore, helpful in understanding the corresponding systematic effects. 

^For a review on twisted mass fermions we refer to e.g. [18] 
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ensemble 

/3 

T/a X (L/a)3 

afii 

a/i^ 

a/i5 

N 

Ns 

Nb 

A30.32 

1.90 

64 X 32^ 

0.0030 

0.150 

0.190 

1367 

24 

5 

A40.24 

1.90 

48 X 243 

0.0040 

0.150 

0.190 

2630 

32 

10 

A40.32 

1.90 

64 X 32^ 

0.0040 

0.150 

0.190 

863 

24 

4 

A60.24 

1.90 

48 X 243 

0.0060 

0.150 

0.190 

1251 

32 

5 

A80.24 

1.90 

48 X 243 

0.0080 

0.150 

0.190 

2449 

32 

10 

A100.24 

1.90 

48 X 243 

0.0100 

0.150 

0.190 

2493 

32 

10 

A80.24S 

1.90 

48 X 243 

0.0080 

0.150 

0.197 

2517 

32 

10 

A100.24S 

1.90 

48 X 243 

0.0100 

0.150 

0.197 

2312 

32 

10 

.625.32 

1.95 

64 X 32^ 

0.0025 

0.135 

0.170 

1484 

24 

5 

.635.32 

1.95 

64 X 32^ 

0.0035 

0.135 

0.170 

1251 

24 

5 

.655.32 

1.95 

64 X 32^ 

0.0055 

0.135 

0.170 

1545 

24 

5 

575.32 

1.95 

64 X 32^ 

0.0075 

0.135 

0.170 

922 

24 

4 

585.24 

1.95 

48 X 243 

0.0085 

0.135 

0.170 

573 

32 

2 

515.48 

2.10 

96 X 48^ 

0.0015 

0.120 

0.1385 

1045 

24 

10 

530.48 

2.10 

96 X 48^ 

0.0030 

0.120 

0.1385 

474 

24 

3 

545.32SC 

2.10 

64 X 32^ 

0.0045 

0.0937 

0.1077 

1887 

24 

10 


Table 1. The dynamical Nf = 2 + 1 + 1 simulations that are included in our investigations. The 
notation that is used to label the ensembles is the same as in [15]. In addition to the value of /3, 
the lattice volumes and the bare quark masses fii, fia-, we give the number of configurations 
N, the number of stochastic samples Ng and the bootstrap block length Nb that were used in the 
determination of flavor singlet quantities [26-28]. 


P 

ro/a 

a [fm] 

Zp/Zs (Ml) 

Zp/Zs (M2) 

1.90 

5.31(8) 

0.0885(36) 

0.699(13) 

0.651(6) 

1.95 

5.77(6) 

0.0815(3) 

0.697(7) 

0.666(4) 

2.10 

7.60(8) 

0.0619(18) 

0.740(5) 

0.727(3) 


Table 2. Chirally extrapolated values of r^/a and Zp/Zg from the two methods Ml and M2 as 
discussed in [25] for our dynamical simulations. 


2.2 Quenched Action 

For the determination of Xoo hr the pure YM theory, we have generated four quenched 
ensembles at four different values of the lattice spacing. For the gauge action, we used 
again the Iwasaki action - Eq. (2.5). 

We emphasize that, in order to test the Witten-Veneziano formula, we tried to achieve 
a very similar setup for the quenched and dynamical situation - we took the same gauge 
action, we matched the physical volume and took a fixed value of rg/i in the valence Dirac 
operator used for spectral projectors, equal to VQfj. of the dynamical simulations^. The 

®We do not have available the value of the renormalization constant Zp for the quenched simulations, 
so this implies only approximate matching of valence quark masses, assuming that with the same gauge 
action Zp is not very different. 


- 5 - 


















/3 

T/a X (L/a)3 

ro/a 

a [fm] 

o/x 

roH 

hue 

2.37 

40 X 20^ 

3.59(2)(3) 

0.1393(14) 

0.0087 

0.0312(3) 

0.158738 

2.48 

48 X 243 

4.28(1)(5) 

0.1182(14) 

0.0073 

0.0309(4) 

0.154928 

2.67 

64 X 32^ 

5.69(2)(3) 

0.0879(6) 

0.0055 

0.0314(2) 

0.150269 

2.85 

80 X 40^ 

7.29(7)(1) 

0.0686(7) 

0.0043 

0.0313(3) 

0.147180 


Table 3. All relevant parameters of the pure gauge ensembles for f3 = 2.37, 2.48, 2.67 and 2.85. 
The errors of ro/a correspond to statistical and systematic uncertainties, respectively. 


details of determining Kc needed for 0(a) improvement can be found in A. All the details 
of the quenched simulations are summarized in Tab. 3. 

The quenched simulations have been performed using the HMC algorithm implemented 
in the tmLQCD package [29]. The usage of this algorithm might introduce somewhat 
larger autocorrelation times compared to e.g. the heatbath with overrelaxation algorithm. 
However, compared to the dynamical simulations, the generation of quenched ensembles 
was still only a small effort. 
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3 Topological susceptibility in the pure Yang-Mills theory 
3.1 The spectral projectors method 

In order to compute the topological susceptibility {xtop) ia the pure YM theory (for this 
case, we will denote it specihcally by Xoo), we used the method of spectral projectors [14], 
This method was originally applied to the pure gauge theory in Ref. [30] and then to the 
dynamical case for twisted mass fermions in Ref. [31]. 

We introduce the dehnition of the topological susceptibility in terms of the spectral 
projector Mm and refer to the original papers for further details about the method [14, 30]: 

Xtop = ^ (Tr{ 75 lRM}Tr{ 75 M|f}) , (3.1) 

where the ratio of renormalization constants Zp/Zs differs from unity due to the use of 
non-chirally symmetric Wilson-type fermions. The spectral projector Mm is an orthogonal 
projector to the subspace of fermion helds spanned by the eigenvectors of the massive 
Hermitian Dirac operator D^D (with D = Dw + in our case) with eigenvalues 

not larger than M^. To achieve O(a^) scaling towards the continuum limit, the renormalized 
value of the threshold M, denoted by Mp, has to be hxed for all ensembles. The value 
of Mp can be chosen arbitrarily, but it is advisable to avoid too small values (close to the 
renormalized quark mass) and choose aMp <C 1 to avoid enhanced cut-off effects [14]. 

In practice, we compute the following spectral observables: 

1 ^ 

' k=l 
1 ^ 

^ k=l 

1 ^ 

c = ^ '^{^MVk, Ib^MVk) , 

' k=l 

where N is the number of stochastic sources rjk {k = 1,..., N) used for the construction 
of the spectral projector, solving the Dirac equation D -|- = rjk an appropriate 

number of times. The above observables are directly related to the right hand side of 
Eq. (3.1) through the following expression: 

(Tr{75Ml^}Tr{75Ml,}) = (C^) - ^ = . (3.5) 

Since the use of two independent sets of sources to compute the square of C would duplicate 
the cost of the calculation, we need to introduce the correction given by {B) /N to correct 
for the bias introduced by the usage of the same set of stochastic sources to compute the 
square of C. 

In order to compute the ratio of renormalization constants ZpjZp^ we also use the 
spectral projector method, since the low cost of the calculation allows us to obtain reliable 


(3.2) 

(3.3) 

(3.4) 
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estimates of ZpjZs not available a priori, since the quenched ensembles were generated 
specifically for this project. The ratio can be obtained through the following expression: 






(3.6) 


This application was first proposed in Ref. [14] and applied for dynamical simulations of 
twisted mass fermions in Ref. [32]. In the latter reference, a comparison of the result given 
by spectral projectors and more standard methods like RI-MOM [33] and x-space [34] was 
done and compatible results were found. A study of discretization effects and finite volume 
effects was also presented. 


3.2 Computation of ZpjZs 

As we mentioned in previous section, in order to compute the topological susceptibility 
using twisted mass fermions, the final value has to be renormalized using the ratio ZpjZs^ 
which is only equal to unity in a chirally symmetric theory. 


/3 

ZpjZs 

2.37 

0.680(1)(27) 

2.48 

0.707(1)(19) 

2.67 

0.752(1)(7) 

2.85 

0.787(1)(3) 


Table 4. Results for ZpjZs values using spectral projectors for the quenched ensembles. The 
errors quoted correspond to, first, statistical and, second, systematic uncertainty (coming from the 
dependence on M). 


In Fig. 1, the value of ZpjZs for different values of tqM for the four quenched ensembles 
listed in Tab. 3 is shown. We follow the strategy presented in Ref. [31] to extract the value 
of ZpjZs for each ensemble, i.e. we take the central value to be the one corresponding to 
some large value of vqM, chosen as 1.5, in order to avoid contamination by non-perturbative 
effects appearing for small tqM. However, due to cut-off effects, even at large tqM, we 
observe some residual dependence on vqM. To quantify the systematic error coming from 
this dependence, we consider the spread of ZpjZs in the interval vqM G [1,2], i.e. this 
systematic uncertainty is taken as the larger of the deviations of ZpjZsiroM = 1.5) with 
respect to ZpjZs{roM = 1) and ZpjZs{rQM = 2). As expected, this uncertainty drops 
down significantly with decreasing lattice spacing. Thus, we obtain the results presented in 
Tab. 4. The first error corresponds to the statistical error given by the spectral projectors 
method, where autocorrelations were taken into account. The second error corresponds 
to the systematic error introduced in our calculation through the residual dependence of 
ZpjZs on M. 

3.3 Continuum limit of Xoo 

The main aim of this work is to test the Witten-Veneziano formula. To this end, we need 
to compute a reliable continuum limit of the topological susceptibility in the quenched 
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I-1-1 
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TqM 
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Figure 1. Results of Zp/Zs as a function of tqM for quenched ensembles. The straight lines 
correspond to the final result and the shaded areas to the systematic error. For further details we 
refer to the text. 


case. The physical conditions, such as the volume or the action used, were matched to the 
situation of the dynamical simulations used to compute the masses. 

Following this strategy, we generated the first ensemble using the Iwasaki gauge action 
at /3 = 2.67 with a value of ro/a matching the one of the dynamical ensemble at (3 = 1.95. 
The other characteristics, such as the volume was matched to the ensemble B55.32 (see 
table 1) with a 32^ x 64 lattice and the value of the valence quark mass for the spectral 
projector method was set to afj, = 0.0055. 

To obtain the continuum limit we needed to extend the simulations to different lattice 
spacings. Thus, we kept the physical volume constant by imposing a constant the physical 
extent of the lattice. In particular, if we assume tq = 0.5 fm, the physical volume was kept 
at L ~ 2.8 fm. Similarly, for the quark mass, we demanded the product to remain 
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p 

(^) 

Ant 

{B) 

Ant 

(C) 

Ant 

2.37 

79.4(2) 

0.5(1) 

17.51(5) 

0.5(1) 

0.19(19) 

0.5(1) 

2.48 

78.7(2) 

0.5(1) 

22.00(8) 

0.5(1) 

0.29(32) 

0.4(1) 

2.67 

78.5(3) 

0.5(1) 

29.4(2) 

0.5(1) 

-0.61(64) 

0.7(2) 

2.85 

78.1(4) 

0.5(1) 

36.5(2) 

0.4(1) 

0.93(93) 

0.8(3) 


Table 5. Results of {A) , {B) and (C) their corresponding values of Tint for quenched ensembles 
at four different lattice spacings. Errors are statistical only. Tint is in the units of the taken step 
between measured configurations, which was increasing for increasing /? (see also Tab. 6). 

invariant. As already pointed out in the previous section, we do not have any estimate of 
Zp and hence we can only match the bare product assuming that the changes in Zp 
are small. All the exact values for each ensemble can be found in Tab. 3. 

As we mentioned in the introduction, the spectral projector method requires a mass 
input parameter M. In principle, the continuum limit result should be independent of this 
value, as long as it is kept fixed for all the ensembles entering the calculation. As the value 
of the renormalization constant Zp needed to renormalize the input parameter M is not 
available, an alternative strategy was followed to match the values of Mp. We used the fact 
that the mode number remains constant at a fixed value of Mp and in a constant physical 
volume [14], i.e. 

oi _ / vin2 

02 \h'2ni 

where o* is the lattice spacing and rii the number of lattice points of the ensembles. Hence, 
the condition In Tab. 5, the values of the mode number given by the spectral observable A 
are shown. 

Once we know the value of the ratio ZpjZg for each /3 value, we are prepared to 
compute the topological susceptibility Xtop- 

The continuum limit is a crucial aspect of lattice QCD. In fact, the rate at which a 
quantity approaches the continuum limit plays a fundamental role in lattice calculations, 
since it is directly related to the accuracy of the final computation of physical observables. 

The twisted mass action at maximal twist guarantees the 0(a) improvement for on 
shell observables [20]. The problem can arise if our observable is affected by short distance 
singularities, since it can spoil the, otherwise guaranteed, O(a^) scaling. 

Our observable, as discussed in detail in Refs. [35, 36], is affected by short distance 
singularities. However, all terms linear in the lattice spacing that are consequently arising 
in the Symanzik expansion vanish at maximal twist [35]. Thus, our observable remains 
0(a) improved, guaranteeing an 0{a^) scaling towards the continuum limit. 

Consequently, we have performed a linear extrapolation in (a/ro)^ to extract the con¬ 
tinuum limit which leads to the following continuum result: 

roXoo = 0.049(6)stat+sys, (3.8) 

where the error is dominated by statistical uncertainty, but takes into account also the 
systematic ones (combined in quadrature with the statistical errors). In Fig. 2, the results 


(3.7) 


10 










{a/rof 


Figure 2. Continuum limit extrapolation of Xoo 


as a function of (a/rg)^ for our quenched ensembles. 


/? 

Ns 

Atraj 

^meas 

r> 

J’oTtop 

2.37 

6 

78000 

769 

25.9(1.3) 

0.0250(14)(22)(8) 

2.48 

6 

98000 

412 

32.8(2.4) 

0.0295(24)(16)(8) 

2.67 

6 

132000 

157 

47.5(5.1) 

0.0376(45)(7)(8) 

2.85 

6 

380000 

119 

60.9(7.9) 

0.0485(70)(4)(18) 


Table 6. Results of Xtop for the continuum limit in the pure gauge theory. Ns represents the 
number of stochastic sources whereas fVmeas is the number of evaluated independent configurations 
and iVtraj gives the corresponding number of thermalized MC trajectories. The errors quoted for 
Xtop are statistical, coming from Zp/Zs and from rg/a, respectively. 


of the topological susceptibility for different values of the lattice spacing together with 
the continuum limit extrapolation are plotted. All the intermediate and hnal results are 
compiled in Tab. 5 and Tab. 6, respectively, as well as other relevant details. 

The continuum value for the topological susceptibility can be compared to earlier results 
using chiral invariant lattice fermions [37], ro4xoo = 0.059(3)stat and from spectral projector 
methods which avoid short distance singularities [30], roi^Xoo = 0.061(6). We note that 
our result is in agreement with the approach used in [38], TqXoo = 0.0524(7)stat(6)sys) which 
has been obtained from a combined continuum and inhnite volume limit. For the older 
result in [37] only statistical errors have been quoted and hnite size effects could not be 


11 ~ 














-30 -20 -10 0 10 20 30 -30 -20 -10 0 10 20 30 

C C 


Figure 3. Histograms of the observable C for four quenched ensembles at /3 = 2.37 (a.), /3 = 2.48 
(b.), 13 = 2.67 (c.) and /? = 2.85 (d.). 

resolved within the statistical accuracy. However, the lattices in the present study exhibit 
a signihcantly larger volume (the physical value of L is about a factor two larger), while 
the range of lattice spacings is similar which might explain the difference in the hnal values 
for TqXoo- In general, older studies [39-42] tend to favor larger values for the topological 
susceptibility as has already been remarked in [38]. 

3.4 Autocorrelations 

We close this section by a discussion of autocorrelation effects in the quenched simulations, 
which might affect signihcantly the errors. 

The topological charge is a quantity highly affected by autocorrelations when the con¬ 
tinuum limit is approached due to appearance of topological barriers, i.e. the transitions 
between topological sectors are suppressed. For this reason, particularly long Monte Carlo 
simulations are needed in order to guarantee that all topological sectors are sampled ade¬ 
quately. 

The topological charge is expected to follow a Gaussian distribution centered at zero in 
the large volume regime. Consequently, due to the fact that the stochastic observable C is 
closely related to the topological charge Q, we expect the same behavior for the distribution 
of C and, equally, (C) = 0. 
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Figure 4. Monte Carlo histories of the observable C for four quenched ensembles at /3 = 2.37 (a.), 
/3 = 2.48 (b.), /3 = 2.67 (c.) and /3 = 2.85 (d.). 

In Fig. 3, the histograms of the stochastic observable C are shown for all the qnenched 
ensembles at fonr different valnes of the lattice spacing. In all cases we were able to construct 
a histogram compatible with a Ganssian distribntion within errors. 

Fig. 4 shows the Monte Carlo history of the observable C again for all the ensembles. It 
is clear, even from visnal inspection of the plots, that antocorrelations become significant for 
the finest lattice spacing. Therefore, even thongh we generated a few times more trajectories 
for this ensemble than for the other ones, we only obtained an estimate of the topological 
snsceptibility with a mnch higher statistical error, since we had to take a step of aronnd 
3200 trajectories between measnrements to obtain Tint compatible with 0.5, i.e. we had 
fewer than 120 independent gange field confignrations (see also Tab. 6). 
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4 Fermionic contributions to the Witten-Veneziano formula from dy¬ 
namical simulations 


Having established a reliable way to extract Xoo from the pure Yang-Mills theory, we still 
need to tackle the l.h.s. of Eq. (1.2). In this section, we discuss how to compute the relevant 
masses and the flavor singlet decay constant parameter /q. The computations in the flavor 
singlet sector are very demanding, because of large contributions from quark disconnected 
diagrams. Moreover, the computation of /o on the lattice turns out to be technically 
involved, due to the fact that the axial vector matrix elements cannot be computed with 
sufficient signal-to-noise ratio. Therefore, we need to consider pseudoscalar matrix elements 
and resort to yPT to relate them to the desired axial vector ones. Further complications 
arise from the fact that the interpolating operators on the lattice are dehned in the so- 
called quark flavor basis, whereas decay constants and the corresponding matrix elements 
are commonly dehned in the singlet-octet havor basis. This issue can be circumvented in 
the framework of yPT as well. 


4.1 Computations in the flavor singlet sector 

The masses for rj and -q' mesons are computed as described in [26, 28]. In particular, we 
employed the method of subtracting excited states in the quark connected part of the cor¬ 
relation functions [43] . Together with the application of the one-end trick for the evaluation 
of quark disconnected diagrams in the light quark sector [44], this yields a substantial im¬ 
provement of the resulting statistical errors [28, 45]. Here, we will only outline some basic 
details that are relevant for the discussion of the flavor singlet decay constant parameter /q 
in the next section. We use pseudoscalar operators in the physical basis: 


j,0,phys ^ ^ 

J)±,PhyS ^ , 


(4.1) 

(4.2) 


where 'ipi, ijji and iph, ijjfi refer to degenerate light and non-degenerate heavy quark doublets, 
respectively. The doublet structure and the flavor projector (1 ± 'r^)/2 is required due 
to the twisted mass formulation and the flavor projector in the second lines allows to 
consider the non-degenerate charm and strange components separately. At maximal twist 
the corresponding operators read 




= O; , 

(4.3) 

^ ' h 1 

(4.4) 


with S the scalar density. In the following, we will neglect the charm operator, as it does 
neither contribute to the q nor the q' within errors, i.e. we consider only the strange com¬ 
ponent ^ pi™-, Pile mixing between the flavor non-singlet scalar and pseudoscalar 

currents in the heavy sector of the twisted basis introduces a relative factor of ZgjZp under 
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renormalization: 


= \zpxh 


(4.5) 


Zs I . i 
r -^75r 

Zp 


Xh ■ 


Similarly, we apply a factor ZsjZp instead of Zs to the light operator allowing to 

pull out a global factor of Z|> from the resulting, renormalized correlation function matrix: 


C{t) = ZlC{t), 


(4.6) 


where 

< 5f’‘'"(t)5f’‘'"(0) > < 5f’*™(t):Pi™(0) > 

< lP*”"(t)5f’*’”(0) > < > 

is build from operators that are renormalized up to Zp and apart from that 

correspond to the ones in Eqs. (4.3,4. 5). 

When building this matrix, we subtract the excited states in the connected contribu¬ 
tions and apply the required factors of Zs/Zp, as given in Tab. 2. We ignore the global 
factor Zp as it will cancel analytically in all observables. Then, we proceed by solving the 
generalized eigenvalue problem [46-48]: 

C (t) (t, to) = (t, to) C (io) (iTo) , to<t, (4.8) 




which gives access to masses through 

(t, to) _ exp -h exp (T - t)) 

A(’^)(t -|- 1, to) exp (—(t -|- 1)) -|- exp (—ml"-) (T — (t -|- 1))) 

and amplitudes 

EC.,(t)uf) (t,to) 

4") = ■■ 7 ■ ■■ - ■ , (4.10) 

J (ul"-) (t, to), C (t) (t, to)) (exp (—m(”)t) ± exp (ml"-) (t — T))) 

respectively. Since masses are renormalization group invariants, they are not affected by the 
actual choice of ZgjZp^ which is only relevant for the computation of amplitudes leading to 
different lattice artifacts for the methods Ml and M2. Note that the amplitudes computed 
in this way are only renormalized up to a factor of Zp as well. 

In Tab. 1, we give the number of configurations iV used on each gauge ensemble to 
extract observables in the flavor singlet sector. Errors are computed using bootstrapping 
with 1000 samples and we have used blocking to deal with autocorrelations. The numbers 
of configurations per block W are given in Tab. 1 and have been chosen to correspond to a 
block length of at least 20 HMC trajectories. In addition, we give the number of stochastic 
samples that has been employed in the computation of quark disconnected diagrams. It 
has been chosen such that the resulting statistical errors are dominated by gauge noise. 

The results for all ensembles are given in Tab. 7. Clearly, the statistical error for a test 
of the Witten-Veneziano formula will be dominated by flavor singlet related quantities as 
can be inferred from the relative errors on the masses. This is caused by large contributions 


- 15 - 






ensemble 

aMps 

aMx 

aMri 

aMrj/ 

A30.32 

0.12358(30) 

0.25150(29) 

0.2800(55) 

0.480(17) 

A40.24 

0.14484(44) 

0.25884(43) 

0.2834(35) 

0.427(14) 

A40.32 

0.14140(30) 

0.25666(23) 

0.2809(28) 

0.458(22) 

A60.24 

0.17277(48) 

0.26695(52) 

0.2908(46) 

0.471(16) 

^80.24 

0.19870(35) 

0.27706(61) 

0.3004(17) 

0.479(23) 

^100.24 

0.22127(32) 

0.28807(34) 

0.3063(24) 

0.454(15) 

A80.24S 

0.19822(33) 

0.25503(33) 

0.2686(30) 

0.463(14) 

A100.24S 

0.22118(33) 

0.26490(74) 

0.2763(13) 

0.518(30) 

R25.32 

0.10685(43) 

0.21240(50) 

0.2348(41) 

0.414(19) 

R35.32 

0.12496(45) 

0.21840(28) 

0.2363(23) 

0.435(24) 

R55.32 

0.15396(31) 

0.22799(34) 

0.2469(28) 

0.480(26) 

R75.32 

0.18036(39) 

0.23753(32) 

0.2537(23) 

0.430(23) 

R85.24 

0.19373(64) 

0.24392(59) 

0.2640(28) 

0.453(28) 

D15A8 

0.06954(26) 

0.16897(85) 

0.1891(64) 

0.295(21) 

D30A8 

0.09801(25) 

0.17760(23) 

0.1969(70) 

0.288(21) 

D45.32SC 

0.11991(37) 

0.17570(84) 

0.1898(22) 

0.276(15) 


Table 7. Results for meson masses on all dynamical ensembles [15, 26, 28]. In addition to the 
relevant masses for the Witten-Veneziano formula in Eq. (1.2), we also give the corresponding values 
of the charged pion mass Mpg. 


from quark disconnected diagrams in the flavor singlet sector which are intrinsically more 
noisy compared to connected contributions and also introduce sizable autocorrelations in 
case of the rj' mass. 

For the computation of the kaon masses, we refer to Ref. [15, 26]. The values are also 
listed in Tab. 7, together with the corresponding charged pion masses Mps. Note that 
some of these values were calculated at smaller statistics than the flavor singlet sector (c.f. 
Tab. 1). Nevertheless, their relative statistical error is still at least one order of magnitude 
smaller than for quantities in the flavor singlet sector. 

4.2 Treatment of /o 

Consider the general dehnition for the decay constant fp, ^ of any pseudoscalar meson P: 

(0|A“(0)|P(p))=z/>^, (4.11) 

where (0) denotes the axial vector current with flavor structure denoted by the index a. 

In the charged sector for mass-degenerate, light quarks, the axial vector current in the 
physical basis transforms into the vector current in the twisted basis at maximal twist. This 
feature can be exploited together with the PCVC relation to derive an expression for the 

^The normalization of decay constants in this work has been chosen consistently s.t. /ps ~ 130 MeV 
holds for the decay constant of the charged pion. 
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charged pion decay constant /pg [49-51]: 


/ps 


2w- 


pr w- 


±\ 


Ks 


(4.12) 


where a = 1,2. Since this formula depends on a pseudoscalar matrix element, it can be 
calculated with very high statistical accuracy. Moreover, this relation allows one to compute 
/ps in the twisted mass formulation without the need for any renormalization factor at all, 
because a factor Zp^ for the bare quark mass /r; cancels the factor for the pseudoscalar 
matrix element in the twisted basis. 

Similar considerations apply for the kaon sector, although a complication arises from 
the fact that one has to employ interpolating operators made of light and heavy quarks. 
The relation for the kaon decay constant fx in the twisted mass formulation is given by: 


Ik 


(W + ^J‘s) 


{0\V, 


,tm 

neutral 


Ml 



(4.13) 


where = \{^{-XdXc + XdXs) + XdilbXc + XdilbXs)- In this case the factor Zs/Zp 

is required again, because it enters fig (c.f. Eq. (2.4)) as well as relative factors due to 
mixing between scalar and pseudoscalar currents. 

Assuming exact isospin symmetry and neglecting possible contributions from the charm 
quark and mixing with further states such as glueballs, the most general parametrization 
for decay constants of the rj, r/' system reads 



/ /8COS(/>8 -/osin(()o 
y /s sin (p8 fo cos (po 


= E{(j)s,4>o)diag{f8, fo) . 


(4.14) 


The choice of Aj), A® together with P = ry, r/' in Eq. (4.11) defines the so-called singlet-octet 
basis. 

Moreover, employing xPT, it is possible to relate the decay constant parameters in 
the r]~r]' system to the remaining octet decay constants /jr, fx and a low energy constant 
Ai = 0{1/Nc) occurring at next-to-leading order in the chiral expansion [11, 12, 52] by: 

/o = (/°) + (/°') = 2 +/.^) + A-i/^ , (4A5) 

fi = ifvf + (4)' = 3 (4/i - fl) , (4.16) 

foh sin(</>8 - </>o) = /°4 + - fn) ■ (4.17) 

In the chosen basis and to the given chiral order, the additional, OZI-violating corrections 

are specific to the singlet sector, i.e. the term ~ Ai = O {^c^) affects neither /§, nor 
any of the angles (po, fs, but only the parameter fo- From the last relation, it can be 
inferred that in the octet singlet basis, the difference between the two angles (po, cp is given 
by SU(3)p’-violating effects, leading to the expectation 


<Po - (p8 
<Po + (p8 


<^1. 


(4.18) 
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The fact that there are different types of contributions to the mixing, i.e. SU (3)^.- 
breaking and OZI-violating effects ~ Ai, can be exploited in order to choose a basis in 
which the two resulting mixing angles do not exhibit a sizable splitting. To this end, one 
introduces the quark flavor basis with the axial vector currents and A® replaced by the 
combinations 


^ (^7m75'« + ^7^75^) , (4.19) 

=\/f =^7 m75S , (4.20) 

in which the light quarks and the strange quark contributions are disentangled. This is the 
reason why this basis is more convenient in lattice simulations, as this is also the preferred 
flavor structure for interpolating operators on the lattice, allowing one to directly access 
the corresponding matrix elements. In exact analogy to the singlet-octet basis, this basis 
again allows for a parametrization in terms of two decay constants and two mixing angles 



= ^{4’i,4>s)diag{fi, fs) , 


(4.21) 


where the mixing matrix H has the same form as the one defined in Eq. (4.14). In this 
basis, the relations between mixing parameters in the rj-r]' system and f-^, fx, Ai read: 

fi = (Af + (40' = 4 + ^Ai/2 , (4.22) 

fs = ( 4 )' + ( 4 )' = 2/i - 4 + ^Ai4 , (4.23) 

fifs sm{(f>i - (j)s) = Affi -k 4 ' 4 = ^Ai/4 (4.24) 


The most important feature of the quark flavor basis becomes manifest in the last expression, 
which is now entirely given by an OZI-violating contribution ~ Ai = O (6), amounting to 
additional suppression for the difference \4>i — (j)s\ compared to {(po — 4>s\, which is given by 
SU(3)ir-breaking effects. Besides, in the SU (3)^;. symmetric case, the angles (pi ~ (pg take 
the value (psTJ{ 3 )p — arctan\/2, and, hence, their numerical value is not expected to be 
small. Therefore, one expects: 

(pi (Ps 

(pi + (ps 

in the quark flavor basis, which has been numerically confirmed in a previous lattice study 
[28]. This feature allows one to consider a simplified mixing scheme in the quark flavor 
basis with only one angle (p 


<C 1 


(4.25) 



E{(p)diag{fi, fs) -kO(Ai) , 


where H {<p) = E (cp, cp). 


(4.26) 
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As mentioned at the beginning of this section, we are restricted in our simulations to 
pseudoscalar operators for practical purposes, because the signal-to-noise ratio for axial 
vector operators turns out too small for a direct computation of the relevant observables. 
However, one may instead consider pseudoscalar matrix elements in order to retrieve infor¬ 
mation on the mixing parameters. This is possible due to the relation between axial vector 
and pseudoscalar matrix elements, which is given non-perturbatively by 

^ (x) 2Mr“i75V> (x) + , (4.27) 

where for a = 0 we have ) and uj (x) denotes the winding number 

density. 

Nonetheless, the anomaly equation of QCD itself is not sufficient for any practical pur¬ 
poses here, mainly because it requires knowledge of an additional matrix element involving 
the topological charge density. Therefore, one needs to gain further insight on how the 
pseudoscalar matrix elements are linked to the mixing parameters. This is again achieved 
by the use of xPT. Consider pseudoscalar currents in the quark flavor basis in analogy to 
Eqs. (4.19),(4.20): 

^ [mj^u + dij^d) , (4.28) 

V2 

P® = siyss , (4.29) 

and the corresponding matrix elements for pseudoscalar mesons P that are given by 

4 = 2mi (0| P* |P) . (4.30) 

To leading order, one can make contact with the quark flavor basis parametrization for axial 
vector matrix elements in Eq. (4.21), i.e. obtain an expression for hp in terms of decay 
constants /;, fs, the mixing angle cj) and octet meson masses [12]: 

Again, formally higher order, OZI-violating contributions are neglected in this expression, 
as demanded by the so-called FKS-scheme, which allows for the determination of process- 
independent mixing parameters [12, 52, 53]. 

Finally, we can combine Eqs. (4.15,4.16) and Eqs. (4.22,4.23) to write down leading 
order relations for the desired singlet decay constant /q in terms of the parameters /ps, /k, 
fl and fs, which we compute on the lattice: 

/o' = -7/6/ls + 2/3/1: + 3/2// , (4.32) 

/o' = +l/3//s - 4/3/i + fl + fs , (4.33) 

/o' = +8/3//s - 16/3/i + 3// . (4.34) 

In general, these relations receive corrections of However, since they were derived 

from continuum yPT, they also differ by lattice artifacts of O(a^), if applied to our lattice 
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data. We will in the following exploit this ambiguity to choose a definition for /o that 
exhibits particularly small systematic effects. Although the use of pseudoscalar matrix 
elements and changing the flavor basis requires to resort to xPT, we point out that this is 
not a serious drawback, as the above expressions are in general of the same order in the 
chiral expansion as the Witten-Veneziano formula in Eq. (1.2). In the following, we will 
refer to the three definitions Eqs. (4.32-4.34) of /o as Dl, D2 and D3, respectively. 

Regarding the computation of decay constants in the charged meson sector, we point 
out that while the values of /ps have been recalculated with the current statistics (they 
were first published in [15]), we used less configurations for fx- However, by far the largest 
contribution to the overall statistical error in our analysis stems again from the flavor singlet 
sector. Finally, we remark that a dedicated study of ?7,7/'-related decay constant parameters 
is currently in preparation [54]. 

4.3 Chiral extrapolations 

Since our lattice simulations employ unphysical quark masses, we need to perform chiral 
extrapolations of our lattice data when computing the l.h.s. of Eq. (1.2). In principle, to 
the given order in xPT, this simply corresponds to a constant fit in (roMps)^. However, we 
have to take into account lattice artifacts which might be different for the three definitions 
Eqs. (4.32-4.34). 

Moreover, the dependence on the choice of ZpjZg is non-trivial, as it affects the relevant 
decay constants (besides /ps) directly through the renormalization of the corresponding 
matrix elements, as well as through the relative renormalization factor that enters the 
quark mass [is in Eq. (2.4), which appears in the definition of fK in Eq. (4.13). 

In Fig. 5 we show results for VqXoo calculated from using our results for the relevant 
meson masses as given in Tab. 7. For the plots a.), c.) and e.) in the left column we used 
ZpjZs values from method Ml and for the right ones b.), d.) and f.) those from M2; c.f. 
Tab. 2. The three rows in Fig. 5 correspond to the three definitions of /o Dl, D2 and D3, 
respectively. 

Clearly, the three definitions of /o show different lattice artifacts and systematic effects 
regarding their ZpjZg (and hence m^) dependence. Besides these systematic effects, their 
relative statistical errors differ as well. Within errors, the first definition Dl shows the 
largest lattice artifacts, as well as the most significant dependence on the choice of ZpjZg- 
Applying a constant fit in (roMpg)^ does not provide a good description of the data in 
this case, as can be seen from panel a.). The fitted value is much lower than for any other 
definition and most points are incompatible with the fitted line. However, the data points 
show a trend towards larger values at smaller lattice spacings, i.e. by fitting only to the 
data at the finest lattice spacing value we obtain = 0.043(12)stat) significantly higher 

than the result of the fit to all data points. 

The data points extracted from definitions D2 and D3 lead to a reasonable agreement 
with a constant extrapolation in (roMps)^, with clearly the best fits stemming from defi¬ 
nition D3. In particular, with definition D3, the extrapolated value is merely independent 
on the choice of Z. 
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method 

MIDI 

M1D2 

M1D3 

M2D1 

M2D2 

M2D3 

roXoo 

0.019(1) 

0.046(3) 

0.058(3) 

0.053(3) 

0.060(3) 

0.036(4) 

xVdof 

4.90 

1.71 

2.60 

2.19 

2.04 

1.95 


Table 8. Results from a constant fit for TqXoo and corresponding x^/dof values (dof = 15) for two 
different sets of values of ZpjZs (Ml,M2) and three definitions for /o (D1,D2,D3); c.f. Eqs. (4.32- 
4.34). Errors are statistical only. 


We have also tried to add a term of O(a^) to our fit. It turns out that apart from the 
most extreme case shown in panel a.) of Fig. 5, the data does not allow to resolve such a 
dependence within errors, i.e. the results for the corresponding coefficient are compatible 
with zero. Regarding the dependence on the strange quark mass, we point out that to the 
given chiral order, the computed Xoo is expected to be a constant function of (roM^y)^ as 
well. Although the data shows a residual dependence on the strange quark mass, no clear 
picture arises with respect to its functional form. 

We have collected numerical results from our constant fits for Xoo in Tab. 8 for the three 
definitions of /o and the two sets of values for ZpjZs. In addition, we give the respective 
X^/dof values. Clearly, the fit MIDI is by far worst, which is expected due to the large 
cutoff effects in this case. 

In order to obtain our final result for Xoo from the dynamic simulations, we apply a 
weight to each fit: 

w = l-2\p-0.5\, (4.35) 

where p denotes the p-values corresponding to the values of x^/dof given in Tab. 8 and 
take the average over all six fits, leading to: 

roXoo = 0.047(3)stat(ll)sys • 

The systematic error has been chosen as the mean absolute deviation from the central value 
and should reflect the uncertainties from residual cutoff and strange quark mass effects. 
Since we included the fit MIDI, which suffers from particularly large lattice artifacts, the 
value of the systematic error should be considered a conservative estimate. 

Another possibility to deal with the residual effects of quark mass dependence and the 
lattice spacing is to include additional, higher order terms in the fit function 

/ {roXoc, {roMxf, (a/rof) = r^Xoc+ciiroM^f+C2iroMKf+C3ia/rof , (4.36) 

where c* with i = 1, 2, 3 denotes the free fit parameters. Note that from the point of view of 
chiral perturbation theory including only a subset of the linear terms would be inconsistent 
with power counting. As can be inferred from table 9 this fit model leads to improved 
X^/dof values. However, most of the additional terms are poorly constrained by the data 
and are close or compatible with zero. In general, this leads to much larger statistical 
errors and all results for TqXoo are compatible within errors. Taking the average of all six 
fits weighted by their respective p-value and statistical errors yields 

r-^Xoo = 0.051 (24)stat, 
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method 

MIDI 

M1D2 

M1D3 

M2D1 

M2D2 

M2D3 

roXoo 

0.062(12) 

0.066(20) 

0.039(31) 

0.049(21) 

0.061(26) 

0.030(38) 

xVdof 

3.45 

1.50 

1.34 

1.91 

1.29 

1.23 


Table 9. Results for TqXoo from the linear fit model in Eq. (4.36) and the corresponding x^/dof 
values (dof = 12) for two different sets of values of Zp/Zs (Ml,M2) and three definitions for /o 
(D1,D2,D3); c.f. Eqs. (4.32-4.34). Errors are statistical only. 


which agrees with the previously computed result from constant hts. Since for the ht 
model in E.q. (4.36) all ht results for r^Xoa are compatible with the averaged result within 
statistical errors, we refrain from quoting an additional systematic uncertainty, as it has 
been done for the constant hts. 

We remark that the inclusion of even more complicated terms seems not feasible, as the 
linear terms in the above model are already rather poorly constrained by the data and close 
or even compatible with zero in many cases. Regarding the data point MIDI it appears 
that the introduction of additional terms does still not lead to a good x^/dof value for the 
corresponding ht, although the resulting value for r^Xoc agrees well with the other results. 
Therefore, we conclude that this data point is rather to be considered a statistical outlier 
than revealing any actual effect related to physics. 
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Figure 5. Results for TqXoo calculated from meson masses and /o as a function of (roMps)^. The 
values of /o employed in panels a.), b.) are obtained from Eq. (4.32), for c.), d.) from Eq. (4.33 
and for e.), f.) from Eq. (4.34). The plots in the left and right column correspond to different 
ZpjZs values from Ml and M2, respectively, as listed in Tab. 2. The chirally extrapolated value is 
obtained from a constant fit in (roMpg)^ and has to be compared with the value computed in the 
pure Yang-Mills theory. 
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Figure 6. Results from the pure Yang-Mills theory (YM) and dynamical simulations. Open and 
closed symbols correspond to the two sets of values for ZpjZs (Ml, M2). All three definitions Dl, 
D2 and D3 for /o are included and represented by circle, triangle and diamond symbols, respectively. 
The solid black line represent the final, p-value weighted average from dynamical simulations and 
its statistical error is indicated by the gray band. In addition, the dotted lines correspond to its 
systematic error; see text. 


5 Discussion 


In Fig. 6, we compare the pure YM topological susceptibility to the left hand side of 
Eq. (1.2) computed in Nj = 2 -|- 1 -|- 1 lattice QCD. All results are given in units of the 
Sommer parameter tq. The weighted mean of the dynamical results and its statistical error 
is indicated by the gray band. The systematic uncertainty of the weighted average is shown 
by the dotted horizontal lines. Apart from the outlier MIDI, which is affected by sizable 
cutoff effects, the results from the dynamical simulations are very close to the one from the 
pure YM theory. Also the agreement of the quenched and the averaged dynamical result 
appears to be good within statistical and systematic uncertainties. 

A complication arises when one attempts to convert the results to physical units. It 
stems from the value of rg in the YM theory: while in = 2 -|-1 -|-1 QCD direct contact to 
physical quantities is natural, in the YM theory such a relation is not obvious. In particular, 
the question arises whether or not the values of rg in the YM theory and in = 2 -|- 1 -|- 1 
QCD are expected to be equal. From a general point of view, we do not see a reason for 
this to be true. 

Therefore, one may use the standard value of tq = 0.5 fm to convert the YM result to 
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physical units and obtain: 


X™ = (185.3(5.6)stat+sysMeV)^ 

For the dynamical simulations, on the other hand, we can use the value tq = 0.474(14) fm 
computed in Ref. [25] to convert to physical units. Taking the aforementioned weighted 
average, we obtain: 

= (193.5(6.2)stat(13.4)sysMeV)^ 

where we have included the error on the physical value of tq in the statistical uncertainty. 
We find rather close agreement between quenched (x^'^) and dynamical (xm°) results 
with deviations of only up to 0(10%). This confirms the validity of the Witten-Veneziano 
formula to the given order and for the assumed value of the quenched tq also in physical 
units. 

6 Summary 

In this paper, we have presented a non-perturbative test of the famous Witten-Veneziano 
formula. This formula relates the large mass value obtained for the rj' meson to the anoma¬ 
lously broken axial U{1) symmetry in QCD. It, therefore, provides important insights for 
our understanding of QCD and the generation of masses. 

We have computed the topological susceptibility in the pure YM theory with dedicated 
quenched lattice simulations using the so-called spectral projector method. In particular, 
by using four values of the lattice spacing, we were able to perform a reliable continuum 
extrapolation and, hence, control this major systematic uncertainty. 

For the first time, we have computed the flavour singlet decay constant /o using Nf = 
2 -|- 1 -|- 1 lattice QCD. Together with the rj, rj' and kaon meson masses determined in 
Ref [28], this allowed us to compute also the l.h.s. of the Witten-Veneziano formula. Again, 
lattice artifacts are controlled by using three values of the lattice spacing. By using a wide 
range of light quark mass values, also the SU(2) chiral extrapolation was performed in a 
controlled way. The strange quark mass dependence was found to be not important for the 
Witten-Veneziano formula. 

The comparison of the pure YM topological susceptibility and its Witten-Veneziano 
counterpart from full QCD leads to agreement within errors. This finding provides clear 
evidence for the hypothesis that the large mass of the rj' meson is due the anomalously 
broken axial ?7(1) symmetry in QCD. 
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A Computation of Kc 



Figure 7. Chiral extrapolation of Kc for the quenched ensemble at /3 = 2.48. The red point 
corresponds to the linearly extrapolated value. 


In order to guarantee the 0(a) improvement, we need to tune k to its critical value 
for the valence quarks. To do so, we followed the strategy introduced in [55]. Thus, we 
computed the of k through the evaluation of rUpcac for different values of the quark mass 
O/U. In particular, we imposed rupcac < O.la/r. In Fig. 7, we show a particular example of 
the chiral behavior of Kc- In all cases, we perform a chiral fit considering only the lowest 
masses a/r < 0.01 since the larger masses deviate from the linear behavior. Notice that the 
data is highly correlated and the fit needed to take this correlation into account. 
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